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We investigate the weakening of elastic materials through randomly distributed circles and cracks 
numerically and compare the results to predictions from homogenization theories. We find a good 
agreement for the case of randomly oriented cracks of equal length in an isotropic plane-strain 
medium for lower crack densities; for higher densities the material is weaker than predicted due to 
precursors of percolation. For a parallel alignment of cracks, where percolation does not occur, we 
analytically predict a power law decay of the effective elastic constants for high crack densities, and 
confirm this result numerically. 
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I. INTRODUCTION 

The appearance of cracks is an effective mechanism for 
a mechanical system under load to release elastic energy 
and to relax towards equilibrium. It is therefore not sur- 
prising that aging processes in a broad class of materials 
can lead to the emergence of microcracks that weaken a 
specimen. They do not necessarily lead to complete fail- 
ure, but their presence alters the elastic properties of the 
system. For many practical applications it is therefore 
highly desirable to develop simple, but still precise pre- 
dictions for the resulting elastic properties of a medium 
that contains defects, cracks or other inhomogeneities. 
Cracked material is just one case in the widely dealt- with 
topic of physical properties of heterogeneous mediai»2i2i4. 
The different physical properties to be described encom- 
pass conductive, transport and also elastic quantities^i 
Often, one starts from a coarse-grained picture and aims 
to find an effective description for the heterogeneous mix- 
ture. Much effort has been put into the calculation of 
effective elastic properties of composed media where the 
constituents have different elastic coefficent o 8 ' 9 ' 10 . 

It turns out that the elastic properties of the system 
depend strongly on the positional and oricntational dis- 
tribution of the inclusions. Even different loading paths 
can lead to a different elastic response of the material un- 
der investigation 1 ^. Therefore, special attention has to be 
paid to the underlying assumptions of the cavity distribu- 
tion. It has been shown that the Hashin-Shtrikman-typ 
bounds obtained for the bulk modulus of two-phase ma- 
terials set important restrictions in terms of phase mod- 
uli and volume fractions^^, and improvement to these 
bounds have to involve considerations of statistical de- 
tails of phase distributions. 

The starting point for various effective medium theo- 
ries is the effect of a single impurity or crack inside an 
otherwise homogeneous medium. Also in the following 
work we will implicitly employ the results of Eshclby- 
concerning the clastic fields around and inside a single el- 



lipsoidal inhomogeneity in an infinitely extended, linearly 
elastic homogeneous solid. In fact, this "dilute" limit for 
a single imperfection already gives an expression for the 
effective elastic constants for vanishingly low defect con- 
centration. Higher concentrations of inhomogeneities can 
be treated in the framework of self-consistent or differen- 
tial effective medium theories. Generally speaking, these 
approaches make use of the idea, that a medium which 
already contains inclusions of another "phase" can be ap- 
proximated as a homogeneous material with different ma- 
terial properties, to which then, step by step, additional 
inhomogeneities arc added to reach a finite concentration 
of them. The underlying assumption, that all effective 
properties depend only on the material constants of the 
pure phases and their volume fractions, is of course only 
an approximation, and the quality of the theoretical pre- 
dictions can hardly be controlled. A careful comparison 
to either experiments or numerical calculations of the 
effective properties is therefore highly recommended to 
judge the quality of the different homogenization meth- 
ods. 

The purpose of the present paper is manifold: First, 
it is intended as a numerical check for the analytical es- 
timates for the effective elastic constants. Obviously, all 
schemes mentioned above are approximative in nature, 
and it is one goal of this paper to shed light on the range 
of applicability of the theoretical models. We use both 
finite difference and finite element methods for the nu- 
merical investigations, and the comparison to earlier re- 
sults serves as benchmark for these approaches. This 
will be done for the important case of spherical inclu- 
sions, since rigorous theoretical statements can be used 
to test the numerical methods. This methodological con- 
firmation is essential for the following tests of homoge- 
nization theories concerning the weakening of materials 
through cracks, which is the second main subject of this 
work. As will be pointed out, the effect of percolation 
plays an important role here, and therefore deviations 
from differential homogenization theories, which predict 
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an exponential weakening of the material, are noticeable 
already for moderate crack densities. In this context, the 
only situation where percolation does not play a role is 
that of parallel cracks. The third important subject of 
this paper is the prediction of effective elastic constants 
in such a geometry, which surprisingly turn out to decay 
here according to a power law behavior, in contrast to 
exponential decays that could be expected from related 
situations^. It must be pointed out that this fully ana- 
lytical prediction becomes accurate in the limit of high 
crack densities, and is therefore complementary to con- 
ventional theories. The predictions are confirmed by the 
same numerical methods that have been justified before. 



II. MODEL VERIFICATION: RANDOM 
DISTRIBUTION OF SPHERICAL HOLES 

The first system under investigation is that of a two- 
dimensional isotropic solid in a plane strain situation 
that contains randomly placed circular holes which are 
allowed to overlap. This system has already been inves- 
tigated numerically by Day et alfii. We use this scenario 
to demonstrate the applicability of our numerical method 
to determine the effective elastic constants. 

For N spherical holes of radius r in the solid phase 
with area A, the true void concentration c is related to 
the void area ratio c = Nwr 2 /A according to the relation 



c = 1 — cxp(— c), 



(1) 



which takes into account that the circles can overlap. 

Starting from the exact expression for a single 
inclusion^, low-density expressions for the effective clastic 
constants can be derived in terms of the two-dimensional 
clastic moduli: 

=£ (2D) - 3£ (2D) c + 0(c 2 ), (2) 
-L Z5) ^ (2D) + (l-3^ ) )c + 0(c 2 ) (3) 

with the elastic constants E^- 2D \ v^ 2D ^> of the solid phase; 
this result is attributed to numerous author o 12 i 13 . We 
use here the explicit annotation 2D to emphasize that 
the elastic constants are those of a two-dimensional plane 
strain material, since some peculiarities in the behavior 
of the effective constants are purely attributed to the di- 
mensionality of representation, as will be elucidated be- 
low. The expressions for conversion between 2D and 3D 
are given in Appendix [X] 

The truncation of the above series after the linear term 
already provides a low density prediction for the effective 
clastic constants. Within this effective medium theory, 
the elastic modulus E vanishes for c = 1/3; however, the 
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0.68, and only then 



E should become strictly zero. This deviation already 
shows that the effective medium theory loses its predic- 
tive power for higher concentrations, underestimating the 
true stiffness of the material. 



Using the above low-density expressions ([2]) and ([3]) , we 
can also derive another approximative model for the elas- 
tic constants in the framework of the differential medium 
theory (see also Appendix [B]) . According to equation 
(IB3I). we start from 



dc 

J (2-D) 

dv\ g 



-3E 



(2D) 



1 

1 - 3f 



(2D) 



dc 1 — c 

and obtain as solution 

E% D \c) =E<? D ){l-c)\ 



eff 3 V 3 



(1 



(4) 
(5) 

(6) 
(7) 



Apparently, this model predicts "percolation" for c = 
1, i.e. if the solid phase disappears completely. It is 
obvious that this model therefore must be invalid for high 
cavity concentrations as well, overestimating the elastic 
constants of the heterogeneous system. 

We note that in both approximative theories the ef- 
fective elastic modulus does not depend on the Poisson 
ratio, a behavior that is known to hold exactly^. 

We use a straightforward finite difference method to 
solve the problem numerically. In a discretized rectangu- 
lar system in the yz plane an "order parameter" <j> is set to 
zero at the grid points which are covered by the circles of 
equal diameter, and <f> = 1 in the remaining solid. Then 
the local elastic modulus is set to E {2D \4>) = 4>E^ 2D \ 
and the elastic equilibrium conditions daij/dxj = are 
solved by relaxation. The system is strained and the av- 
erage stress calculated, from which the effective clastic 
constants can be deduced as follows: For a system that 
is strained in z direction and has periodic boundary con- 
ditions in y direction, the average strain (e yy ) vanishes. 
The average diagonal stress components in the system 
(<7y y ) and (a zz ) are measured for this plane strain sce- 
nario. Then the effective elastic constants are determined 
through 

E (3D) _ (^{(Tyy) + (cT zz ))((cT zz ) - (<Tyy)) 



J eff 



,(3£>) 
'eff 



{{(Jyy) + ((T ZZ ))(e zz ) 

{ a yy) 



{(Jyy) + (a zz ) 7 



(9) 



where the average strain (e zz ) is fixed through the bound- 
ary conditions. We typically used systems sizes of 
2048 x 1024 grid points, with up to 1000 circles with 
a radius of 20 grid points. Further details on the clastic 
solver are presented i n 17 ' 18 . 

The dependence of the effective elastic modulus on the 
concentration as predicted by the theories, see Eq. @ 
and Eq. ([6]) , and as obtained by numerical simulations is 
shown in Fig. [TJ The independence of E e g on the Pois- 
son ratio is clearly visible also in the numerics, where we 
checked this explicitly for v {2D ^ = 1/2 and = -0.41 
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FIG. 1: Effective elastic modulus as function of the void con- 
centration c. The plot shows numerical data for different Pois- 
son ratios, as obtained with the present method, in compari- 
son to numerical results obtained by Day et al. Here we used 
different distributions of cracks, and evaluated the effective 
elastic modulus for the two different Poisson ratios using ex- 
actly the same arrangement of cracks; the independence of 
the Poisson ratio is clearly visible. 



FIG. 2: Poisson ratio as function of the void concentration c. 
For v = 1/3, both, the effective medium theory and the differ- 
ential theory show fairly good agreement with the numerical 
results. The predictions from the effective medium theory are 
shown only up to the percolation point c = 1/3. For negative 
Poisson ratios, the differential theory coincides much better 
with the simulations. 



(corresponding to i/ 3£) ) = 1/3 and i/ 3 -°) = —0.7 respec- 
tively); the latter case of an auxetic material is sometimes 
observed e.g. in foams^, and is here only used as an ex- 
treme case to confirm the independence on the Poisson 
ratio. In fact, wc find that for the same random ar- 
rangement of circular holes the elastic constants match. 
Since we wanted to obtain a reasonable statistical av- 
eraging, we also performed repeated runs with different 
initializations. As we increase the void concentration c, 
one can clearly see that the scattering of the data points 
increases for higher concentrations, since larger clusters 
can form which can become comparable to the (finite) 
system size used in the simulations. Also, the relaxation 
time increases strongly with c, thus results for higher 
concentrations are not shown here. In Ref. [ll|, Day et 
al. performed simulations based on an elastic spring net- 
work formulation for system setups analogous to ours. 
The comparison of our numerical results to the simula- 
tion data of Day et al. are also included in Fig. [T] The 
results for the independent numerical approaches are in 
reasonable agreement. In particular, all sets correctly re- 
produce the exactly known low density limit c — > 0. For 
higher concentrations, we obtain a higher effective elas- 
tic modulus than Day et al., and we believe that this is 
a consequence of the considerably larger systems that we 
used. 

Similarly, the effective Poisson ratio agrees well with 
the differential theory, as can be seen in Fig.[2j especially 
in the case of a negative Poisson ratio. Even at the high- 
est densities that were simulated here, we do not observe 
a noticeable deviation from this homogenization model. 

Finally, we briefly remark that the results depend on 
the dimension of representation. Conversion of the re- 



sults for the differential homogenization theory gives ac- 
cording to Eqs. (|A2[) 



E 



(3D) 
eff 



3£( 3D > (c - l) 3 (c(8u^ - 2)(c 2 - 3c + 3) 



3(1 + ^ (3D) )) / {c(Av 



l)(c 2 -3c + 3) -3) 2 



x(z/ 3D ) + l) 



.(3D) 
eff 



c (Ay 



(3D) 



l)(c 2 -3c + 3)-3^ 3D ) 



c(4i/( 3I) )-1)(c 2 -3c + 3) 



(10) 
(11) 



In particular, the effective three-dimensional elastic 
modulus does not have the property of being independent 
of the Poisson ratio. Furthermore, for negative Poisson 
ratios the effective elastic modulus can first increase if the 
material is "weakened" by spherical holes. A similar be- 
havior was reported for cracks in Ref. 0, and here we see 
that this effect is rather generic and results mainly from 
the definition of the elastic constants. Indeed, this coun- 
terintuitive behavior is obviously an artifact of the three- 
dimensional representation that is already contained in 
the low density expressions and not related to a specific 
homogenization scheme. Already for low concentrations 
we get 



(3D) 
eff 



* (a "'(i- (i -i;y )+3> ^>. 

(12) 

which can start with a positive slope for negative Poisson 
ratios. 
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III. RANDOM DISTRIBUTION OF CRACKS 

In this section we investigate a random arrangement 
of cracks in a solid and compare the prediction for the 
effective clastic constants to numerical simulations. To 
that end, we use the same geometry as in Ref. [tl where 
the normal vectors of the planar cracks are located in 
the yz plane, and they arc infinitely extended in x di- 
rection. Therefore, the system becomes again effectively 
two-dimensional, and we restrict our investigations to a 
plane-strain scenario. In the yz plane, all cracks have 
the same length L; here we assume that the orientation 
is random and all angles 9 appear with the same proba- 
bility; in the notation of Ref. 13 this means for the orien- 
tational order parameter P — (sin 2 9) = 1/2. 

We introduce a crack density parameter 



tt(L/2) 2 N 
A ' 



(13) 



where N is the number of cracks per area A in the 
yz plane. The prediction for the effective (three- 
dimensional) elastic constants in the framework of the 
differential homogenization method is for plane strain ac- 
cording toS. 



E 



(3D) 
eff 

,(313) 



:i 



.(3D)) 



[2,(31?) [/ (3r>)) e a]2(' 1 + ^(31?)) 



,(31?) 



1 



/(31?)) e a 



(14) 
(15) 



which predicts an exponential weakening of the material 
with the density parameter a. In particular, the effective 
medium is still isotropic, since there is no preferred orien- 
tation for the cracks, and therefore the elastic properties 
are still fully described by two elastic constants. 

Interestingly, the two-dimensional representation of 
the above result gives simply 



E 



(2D) 
eff 



E^ exp(-a), 



,(21?) 
eff 



,(21?) 



exp(- 



■a), 
(16) 

so both constants decay according to a simple exponen- 
tial decay to zero. Notice in particular that the effec- 
tive Poisson two-dimensional ratio also tends to zero, in 
contrast to the spherical case discussed before, where it 
approaches 1/3. We also mention that here the effective 
elastic modulus does not depend on the bare Poisson ra- 
tio v^ 2D \ Notice that the above conversion implies also 
that the effect of an increase of stiffness with the crack 
density for negative Poisson ratios, that was discussed in 
Ref. |3, is indeed an artifact of the three-dimensional rep- 
resentation, similar to the spherical example discussed 
above. 

We note that the limit E e g = is only reached for 
a — > oo, which means that this theory does not pre- 
dict percolation. However, in reality percolation occurs 
for— a w 4.49, and then a network of cracks penetrates 
the whole system, thus the true effective modulus van- 
ishes. Therefore the differential homogenization method 



0.8 



0.6 



'% 0.4 



0.2 



PL. Theory 


v< 3D >=1/3' 


Numerics 


v< 3D >=1/3 x 


xV 




^ N 


x — 




X X 



0.5 



1 1.5 

a 



2.5 



FIG. 3: Effective elastic modulus as function of the crack den- 
sity for plane strain loading, v^ aD ^ = 1/3, for several random 
distributions (P = 1/2) of cracks of equal length. 
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FIG. 4: Effective elastic modulus as function of the crack den- 
sity for plane strain loading, u^ iD ^ = —0.7, for several random 
distributions (P = 1/2) of cracks of equal length. The initial 
stiffness increase predicted by the differential homogenization 
theory is clearly visible. For higher crack densities, the theory 
overestimates the effective elastic modulus significantly. 



overestimates the true elastic modulus for higher crack 
densities. 

To check the quality of the above analytical predic- 
tions, we investigated the case of randomly oriented 
cracks also numerically for plane strain using finite differ- 
ence relaxation methods, see Figs. [3] and 2J For low crack 
densities, the numerical results agree with the prediction 
Eq. (|14[) , but for higher values they are indeed systemati- 
cally lower due to prospective percolation. Nevetheless, 
the analytical prediction from Ref. [3" can be considered 
as a very good approximation at least for crack densities 
a < 1. 

We also see good agreement for the Poisson ratio in 
this range of a, see Fig. [5] For a negative bare Poisson 
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FIG. 5: Effective Poisson ratio as function of the crack density 
for plane strain loading for several random distributions (P = 
1/2) of cracks of equal length. 
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FIG. 6: Random arrangement of parallel cracks. The aver- 
age crack length is L, the average vertical distance between 
neighboring cracks h. 



ratio, here i/ 3 -°) = —0.7, the numerical results seem to 
indicate that it approaches even a positive value instead 
of just decaying to zero. 



IV. ASYMPTOTIC BEHAVIOR OF PARALLEL 
CRACKS 

From a more general point of view, all setups with ran- 
dom crack orientations have a finite percolation thresh- 
old, even if the probability distribution for the choice of 
the angle is not uniform. The only exception is the case 
that all cracks are parallel; then percolation does not oc- 
cur. Thus only here a nontrivial asymptotic behavior 
exists for high crack densities. It turns out that for this 
special case analytical predictions for the effective elastic 
constants can be made, which become accurate in the 
limit a — > oo, and in this respect they differ fundamen- 
tally from conventional homogenization theories. 

First, it should be noted that in this case the material 
that is pierced by cracks becomes anisotropic, and we 
therefore characterize its elastic properties by the tensor 
Cijkitki- If we assume that in the yz 



c tfkn with a H 



plane all cracks are aligned in y direction (see Fig. [S]), 
it is immediately clear that e.g. c y & yy = c yyyy , since a 
pure stretching in y direction does not open the cracks; 
hence the strain tensor is homogeneous in the material 
and unaffected by the cracks. 

For low crack densities, the effective elastic constants 
were calculated in Ref. [^, and in particular we get 



C3333 = (1 - + 2aP]D- 1 E^ 



(17) 



with 



D = [4a 2 P(l -P)(l 



v 



(3D)\2 



) 2 + 2(l-v^) 2 a + 



+l- 2l / 3D )](l + z/ 3D) ) 



and P = for the parallel arrangement. 

We start with looking at high crack densities, a — > 00: 
two different lengthscales are important for a complete 
description of the problem at hand, the length L of the 
cracks and the average vertical distance h between them. 
For high crack densities a, the vertical distance h be- 
tween neighboring cracks is smaller than the average 
crack length L, and the relation between the two char- 
acteristic length scales can be given through a only, so 
we obtain h ~ L/a. If the cracked body is subjected 
to tensile loading perpendicular to the cracks, the solid 
regions between two cracks can be understood as a thin 
bent plate of a width proportional to L and thickness h. 
The opening of the cracks is the displacement u z . The 
stress of a thin bent plate scales a o 21 i 22 



Eh 3 



1 



(18) 



With this equation, it follows readily that the average 
stress and the opening u z have to scale like 



(u z ) ~ (a zz ) 



(1 - v 2 )L 4 
Eh? 



(19) 



The total displacement is distributed among the opening 
of all cracks, which relax the material around them. Since 
for this loading all other average strain components are 
small—, the average strain (e zz ) is simply given by 



(20) 



Plugging this into Eq. (TIT))) , we finally obtain for the case 
a > 1 



E 



(1 - v 2 )a 4 



(21) 
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In other words, the relevant elastic constant 



decays by a power law, 

Csfas ~ c 3333«" 4 - (23) 

We note that this scaling behavior holds also for situ- 
ations where the cracks can have unequal lengths, dis- 
tributed around the mean value L; details of the distri- 
bution function can affect only the numerical prefactor 
of the above prediction in the limit a — > oo. In addi- 
tion, we also performed simulations for regular arrays of 
cracks. Also, we checked numerically that the scaling be- 
havior holds for random parallel arrangements of cracks; 
the results can be seen in Fig. [JJ This graph shows the 
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a 

FIG. 7: Scaling behavior of the effective elastic constant Cg^ 33 
as a function of the crack density a for a parallel arrange- 
ment of cracks in logarithmic representation. For a regular 
arrangement of cracks, the agreement of the numerical simu- 
lations with thin plate theory is excellent. If the cracks are 
placed at random positions, they still exhibit the same power 
law scaling behavior. 

results for the low density theory, the asymptotic behav- 
ior and numerical simulation data from both finite dif- 
ference and finite element methods^. We used different 
arrangements of cracks to illustrate the scaling behavior: 
First, we took a regular arrangement of cracks, where we 
can rigorously calculate the effective elastic constants for 
a — > oo; this is shown in Appendix [C] Due to the spa- 
tial periodicity it is sufficient to consider a system with 
only a few cracks. We clearly see that both finite differ- 
ence and finite element calculations give the same result. 
The finite element method is computationally more effi- 
cient than the simple relaxation solver; however, the geo- 
metrical description is easier with finite differences, since 
e.g. intersections with boundaries (or overlaps of cracks 
for the random orientation discussed in the pre- 

ceding section) do not require a separate treatment. To 



get clear predictions for the scaling behavior as function 
of the crack density a, we randomly place the cracks in 
the system and solve the elastic problem by finite element 
methods. Then we change the value of a by rescaling the 
height of the system, which means that the arrangement 
of cracks is the same for all points on one curve. The 
correct scaling behavior is demonstrated here for a rela- 
tively small system with only N = 20 cracks. Obviously, 
the specific results depend then on the configuration, and 
only for N — > oo these discrepancies between different ar- 
rangements would disappear. However, the results show, 
that the scaling holds for each configuration (shown here 
for two cases), and therefore it must be correct also for 
the true ensemble average in an infinitely large system. 

The results, in particular the finite difference data for 
small a show the crossover between the low density pre- 
diction (|17p and the asymptotic behavior ([23]). For the 
latter, the numerical prefactor was chosen such that it 
matches the particular case of regular cracks (g = 1/2), 
as explained in Appendix [Cl 



V. SUMMARY AND CONCLUSION 

We investigated numerically the effective elastic con- 
stants for isotropic plane strain media with spherical 
holes, randomly oriented and parallel cracks. In all cases 
we find a good agreement with predictions from different 
homogenization theories, with a better performance of 
differential media theories. The results show clear devia- 
tions from the approximative theories, which are strictly 
valid only for low inclusion densities, since they do not 
correctly account for effects which go beyond mean-field 
approximations. In particular, all discussed models do 
not correctly take into account percolation, which should 
lead to a sharp drop of the effective elastic modulus. The 
only case where percolation does not occur is that of par- 
allel cracks. By scaling arguments we derived analytically 
the scaling behavior of effective elastic constants in the 
limit a — > oo and obtain a power law decay with the crack 
density. This new prediction was confirmed numerically 
using finite-difference and finite-element methods. We 
note that this prediction is complementary to conven- 
tional homogenization theories, as it becomes accurate 
for increasing crack densities. Even though the effective 
clastic constants are already low in this regime, the ob- 
tained results are therefore of principal interest and raise 
the question whether explicit solutions for other situa- 
tions with high inclusion density are also possible. 
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APPENDIX A: CONVERSION BETWEEN TWO- 
AND THREE-DIMENSIONAL 
REPRESENTATION 



As already mentioned above, the dimensionality can 
play a role for the effective elastic constants. We can 
convert the elastic constants of a two-dimensional setup 
to an equivalent three-dimensional plane strain situation. 
The defining equation is Hooke's law, 



1 r . 



HjCTkkl 



(Al) 



which holds for both 3D and 2D; the difference is that in 
the first case all indices run over x, y, z, in the second only 
over y,z. In a plane strain 3D configuration, e xx = 0, we 
have a xx = v(a XJ y + a zz ), whereas this stress component 
does not appear in 2D. Hence the conversion rules for the 
elastic constants are given by 



£(3 



D) 



2D) 



1 + 2t/ 2J3 > 

(1 + Z,(2£>))2 ' 



A 2D ) 



which follows directly from Hooke's law (jAlj) . 



(A2) 





FIG. 8: Top: Sketch of the regular array of cracks that is used 
both for analytical calculations and numerics. The dashed 
rectangle is the "periodic unit cell" in which the elastic prob- 
lem is solved numerically. The dark box visualizes the plate 
that is bent under the applied load, which is shown in the 
lower panel. The deformation of the neutral fiber is denoted 
by z(y). 



APPENDIX B: DIFFERENTIAL 
HOMOGENIZATION METHOD 

Let a system of dimensionless "volume" Vq = 1 contain 
inclusions of a second phase, characterized by the initial 
concentration (volume fraction) cq , which in turn means 
that the concentration of the first phase is 1 — cq. Now, 
a volume dco of the second phase to the original volume 
Vo = 1 is added, leading to a total volume of V = 1 +dc . 
The total volume of phase two has increased to cq + (Icq, 
resulting in a total volume fraction of 

c = C ° + = co + (1 - co)dc + C(dcg). (Bl) 
t + clcq 

The change of concentration of the second phase is there- 
fore dc = (1 — c )dco. Let M e g denote a complete set 
of effective elastic constants (or other quantities of inter- 
est). In the framework of the homogenization methods 
used here, this set should depend on the properties of the 
pure phases and the concentration, M e g = F(Mi, M 2 , c) 
with a universal function F and the obvious relation 



dco to the already homogenized medium with properties 
M e g. Hence we obtain 

M eff + dM eff = F(M 1 ,M 2 ,c + dc) = F{M eff ,M 2l dc ) 



M, 



dF(M 1 ,M 2 ,c) 



dc 



dc 



c=0,M 1 =M eff (c) 



1 



since in the second step the change of concentration is 
dco/(l + dc ) = dc + O(dcl); in the last step the relation 
(|B2[) was used. Here it is important to note that after the 
differentiation first c has to be set to zero, and only then 
Ml = M e ff(c) to be inserted. Therefore, we immediately 
obtain the fundamental equation 



dM, 



eff 



I dF(M 1 ,M 2 ,c) 



dc 



1 



dc 



(B3) 



c=0,Mi=A%(e) 



For slit-like cracks, the volume fraction is zero, and there- 
fore the prefactor (1 — c)~ l disappears and c is replaced 
by the density parameter a. 



M eff = F(M eff ,M 2 ,c = 0). 



(B2) 



In the framework of the differential homogenization 
method the increase of the amount of the new phase from 
c to c + dc is interpreted as the addition of the amount 



APPENDIX C: REGULAR ARRAY OF CRACKS 

To make the preceding scaling arguments in section UVl 
more explicit, we discuss here a regular arrangement of 



8 



cracks, as depicted in Fig. [8] and solve this problem ex- 
actly in the limit a — > oo. The idea is that the displace- 
ment, which is applied to the sample is mainly stored in 
the opening of the cracks, and the material in between 
is only slightly stretched. The region between adjacent 
cracks behaves then as a bent plate (see dark region in 
Fig. |5J), which is thin in the limit R ^ h. Wc note 
that for this regular arrangement the plate length R ap- 
pears here as additional parameter, which is related to 
the gap distance x by L = 2R + x; again, L is the crack 
length which is now assumed to be exactly the same for 
all cracks. Therefore, the additional dimensionless pa- 
rameter g = x/R remains in the final solution, whereas 
for an irregular arrangement of cracks it would be deter- 
mined statistically; finally, it enters only into the numer- 
ical prefactor of the effective elastic constants. 

For the given geometry, the area that is occupied by 
a single crack, N = 1, is A = {L + x)h. Therefore, the 
crack density is 



n (l + g/2) 2 R 

Q' — 

2 l+g h 



(CI) 



The bending of the thin plate is described by the equation 
z ""{y) = 0, since the upper and lower surfaces are stress 
free^i. Each plate is displaced by z(R) = (e zz )h, since 
the total displacement is equally distributed among all 
crack openings. Together with the symmetry conditions 
z'(0) = z'{R) = and the reference value z(0) — 0, wc 
obtain for the coefficients of the general solution z{y) = 



ay 3 + by 2 + cy + d the values b = 3(e zz )h/ R 2 and a — 
—2b/3R. The force per unit length in x direction that is 
required to bend the plate by the given amount is given 
by^i 



F = —- 



Eh 3 ,„ _ Eh A ( 



12(l-i/ 2 ) (l-^)i? 3 ' 
and thus the average stress in vertical direction 



(cr zz ) = 



F 



E(e zz ) /7TX4 (l+ 5 /2) g 



x + R 1 — v 2 V 2 / (l+g) 



(C2) 



(C3) 



From Hooke's law for the effective medium, (a z 

c 3333( £ zz) + . . . follows 



-3333 



\-v 2 \2 



The bare elastic constant C3333 is related to the isotropic 
moduli by 



C3333 — 



£7(1-1/) 



(1 + i/)(1-2i/)' 
and hence get get asymptotically for a — > 00 



- eff \-2v /7TX4 (l+g/2) 8 _ A 

C3333 ~ (l-^) 2 V2/ 1 (l+ 5 )5 a 



3333 



(C5) 



(C6) 
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